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Abstract 

In this work, we study the motions in the region around the equilateral Lagrangian 
equilibrium points L 4 and L 5 , in the framework of the Circular Planar Restricted Three- 
Body Problem (hereafter, CPRTBP). We design a semi-analytic approach based on some 
ideas by Garfinkel in [3]: the Hamiltonian is expanded in Poincare-Delaunay coordinates 
and a suitable average is performed. This allows us to construct (quasi) invariant tori that 
are moderately far from the Lagrangian points L4-L5 and approximate wide tadpole orbits. 
This construction provides the tools for studying optimal transfers in the neighborhood of the 
equilateral points, when instantaneous impulses are considered. We show some applications 
of the new averaged Hamiltonian for the Earth-Moon system, applied to the setting-up of 
some transfers which allow to enter in the stability region filled by tadpole orbits. 


1 Introduction 

For the design of simple transfers in Astrodynamics, Hohmann transfers are widely used. They 
consist basically in a maneuvre in a system represented by a Two-Body Problem (2BP) and 
their solutions, consequently given by Keplerian ellipses. Starting from a circular orbit around 
a main body, a transfer to a different (inner or outer) circular orbit can be achieved with just 
two different impulses of properly defined sense and magnitude (see e.g. my 

In cases where 2BP is not a suitable approximation, a similar approach can be considered 
with a different simplified version of the Hamiltonian representing the system. The aim of 
the method lays on the idea of designing a transfer between orbits that are exact solutions of 
an integrable approximation of the studied system, using a set of impulses. Since, in general, 
physical problems have more than one degree of freedom, a suitable construction of a normal 
form approximating the Hamiltonian of the system is needed. Furthermore, many d.o.f. make 
representations in configuration space inadequate, since they do not give precise hints of the 
time evolution of the orbits at glance. Thus, the baseline is the construction of a normalized 
integrable approximation of the model to study, that should provide: i) analytical solutions for 

* Key words and phrases: Restricted three-body problem, normal forms, Hamiltonian perturbation theory, 
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the motions, and ii) suitable surfaces of sections, tools to be used for the analysis of the effects 
of the impulses that conform the trasfer. 

In the last decades, several semi-analytical results have been obtained in order to ensure the 
stability of the motions of some Trojan asteroids, orbiting around the Lagrangian equilibrium 
points L 4 — L 5 of the Sun-Jupiter system. All those works share a same common structure, 
summarized as follows: each approach is based on an explicit algorithm, that can be translated 
on a computer so as to calculate the expansion of a suitable normal form, providing a good 
local approximation of the complete Hamiltonian of the CPRTBP. Such a normal form is used 
to approximate the orbits of some objects and to prove their stability, provided they are close 
enough to the equilateral Lagrangian points. As far as we know, the wider coverage of the 
tadpo 1(0 orbits around L 4 — L 5 is given in [3], that is based on the Kolmogorov normal form. 
Here, we try to improve those results, by revisiting the approach developed in [^. From that 
article we borrow two main ideas: first, we perform the initial expansions of the CPRTBP 
Hamiltonian in a suitable set of Poincare-Delaunay-like coordinates (while polar coordinates 
were used in i). That particular type of canonical variables allows us to clearly distinguish 
a pair of slowly varying coordinates from those quickly changing their values. Therefore, we 
average the Hamiltonian with respect to the angle related to the fast dynamics. This second 
main idea is implemented here, by using the modern Lie series formalism so as to construct an 
integrable normal form; this allows us to approximate also the tadpole orbits going rather far 
from L 4 — L 5 . In fact, as a major novelty with respect to [3] (that is based on purely analytical 
techniques), we have translated our procedure in some codes, producing suitably truncated 
expansions of the normal form and, therefore, explicit numerical results. 

Having the tools described before, we create an algorithm which is able to test the effect¬ 
iveness of different impulses, after choosing an arbitrary starting point. This approach can be 
applied in many different systems. In particular, we focus our work in the area surrounding the 
Lagrangian equilibrium points L 4 — L 5 of the Earth-Moon system. This kind of stability region 
is interesting since it is located close to the Earth, and could provide a natural trapping zone 
for small bodies as astronomical observatories, obsolet spacecrafts or space debris in general. 
At the end of the paper, we characterize the effects of different instantaneous impulses and we 
choose the best candidates for effective transfers in the framework provided by the normal form 
approximating the CPRTBP Hamiltonian. 


2 Explicit construction of the integrable approximation 

2.1 Initial settings 

Let us introduce the standard framework of the CPRTBP, as done, e.g., in [2]. In such a model, 
the motion of the two biggest bodies (hereafter, the primaries) is not influenced by the third 
one, considered massless, so the orbits of the primaries are Keplerian ellipses and the third body 
moves under the gravitational attraction exerted by the other two. Additionally, we assume that 
both those Keplerian orbits are circular and the massless body is coplanar with the primaries. 
Let us define the heliocentric vectors r^o = Xj — xq with j = 1,2, being xq , xi and X 2 the 
position vectors of the biggest primary, the smallest one and the third body, respectively, in an 
inertial frame. As usual, we also define the units of measures in order to set the rotation period 


^For an introduction to tadpole and horseshoe orbits in the CPRTBP model, see, e.g., § 3.9 of [6]. 
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r = 27r and the gravitational constant Q = 1, and we denote with p and 1 — /r the masses of the 
smallest primary and the largest one, respectively. These settings imply that the semi-major 
axis of the ellipse described by the orbit of rio is equal to 1, while p E (0,1/2). 

For our purposes, it is convenient to adopt the Hamiltonian formalism in the non-inertial 
synodic frame that co-rotates with the primaries. We define the axes so that the largest primary 
is always located at the origin of the synodic frame, while the fixed position of the smallest 
primary is such that rio = (—1,0). Let us first introduce the action-angle canonical coordinates 
(G, r. A, i), which can be seen as modified Delaunay variables for the planar Keplerian problem. 
They are given by 


G = Vo(l - e2) , x = M + g-M' -g' , 

r = Va (l - Vl - , £ = M , 


( 1 ) 


where a, e, M and g are the semi-major axis, the eccentricity, the mean anomaly and the 
longitude of the perihelion of the massless body, being the angle g measured in the inertial frame. 
Moreover, M' and g' denote the mean anomaly and the perihelion longitude of the smallest 
primary, respectively. Thus, A corresponds to the synodic mean longitude. The Hamiltonian 
ruling the motion of the third body in the non-inertial synodical frame can now be written as 
follows: 

n{G, r. A, ^ - G - pF{G, r, a, i) , ( 2 ) 

where the so-called disturbing function pF is such that 


F = 




+ 


rio • ^20 

l|rioP 


(3) 


In order to remove the singularity of the Delaunay-like variables when e = 0 (i.e., for circular 
orbits), it is convenient to introduce canonical coordinates (p, ^,A, ry), similar to the Poincare 
coordinates for the planar Keplerian problem. Thus, let us define 


p — G — 1, A — A, 

^ cos £ , rj = sin £ , 


(4) 


where the values of p are significantly small in a region surrounding the Lagrangian points, for 
instance, in the case of tadpole or horseshoe orbits. In fact, in those cases, it holds true since 
a ~ 1 and e > 0. Let us recall that those variables have been adopted in [2] to study the 
CPRTBP. 

Before constructing a normal form, the starting Hamiltonian ([2]) (in particular, the disturbing 
function ([3])) must be expanded in Poincare-Delaunay-like coordinates {p,^,X,g). Such a non 
trivial operation is described in m, which also includes all the detailed Mathematica codes 
explicitly producing those expansions. This approach ensures that the starting Hamiltonian T-L 
can be written in the following form as a function of the Poincare-Delaunay-like coordinates: 


iL(°’°)(/9,^,A,ry) = '^ zj°\p, + g'^)/2) X,g) , (5) 

l>0 s>l l>0 

where E Tip and £ Fi,sK V / > 0, s > 1, being K a fixed positive integer. VpsK is 

the set of functions such that 








4 


R.I. Paez, U. Locatelli 


• a function g G Vi^sK if the generic terms appearing in its Taylor-Fourier expansion, which 
are of type cos(feA) or sin(A:A), satisfy the fol¬ 

lowing relations about their coefficients: 

Cmi,m2,m3,k = dm3,m2,m3,k = 0 ^hen 2 mi m2 -h m3 / / or \k\ > sK . 

In principle, the expansion ([5|) can be seen as a reorganization of the Taylor-Fourier series giving 
the disturbing function; this is made in a suitable way that allows us to successfully perform the 
construction of the normal form. Furthermore, the criterion for the choice of K is such that the 
size (in any common functional norm) of ... is approximately the same for any 

value of the index 1. In other words, the disturbing function is splitted in terms 0{g), 0{g?‘), 
... by using the Fourier decay of the coefficients for increasing values of the harmonic \k\; let us 
recall that in the present work g is regarded as a fixed small parameter of the system. 


2.2 Averaging the Hamiltonian over the fast angle 

Let us emphasize that the expansion © contains also a non-trivial information about the 
Keplerian part (that corresponds to the whole Hamiltonian if /r = 0). In fact, since = 
{Pi (C^ + ^^)/2) and zf''^ € Vi^ , then one can deduce that zf''^ = 0 when the index I is odd. 
This is in agreement with the expansion of the starting Hamiltonian T-L, when the disturbing 
function is neglected and the actions G and F are substituted according to formula (H]). If we 
explicitely write the hrst main terms of the Keplerian part 


7 ( 0 ) I y{ 0 ) 


(0)_ 3 + 3 

+ ^4 - 2^ 2 2 


P + 


-k rf 


= -| + r-|(p + r)^ 


( 6 ) 


we conclude that the angular velocities have different order of magnitude 

; dU 


A = — ~ 0 , 
op 


dr 


~ 1 


(7) 


because /r <C 1 and the values of the actions p and F are small in a region surrounding the 
Lagrangian points. Thus, A can be seen as a slow angle and £ as a fast angle. Then, this 
motivates to average the Hamiltonian over the fast angle (see, e.g., § 52 of lU), in order to 
focus mainly on the secular evolution of the system. We remove all the terms depending on 
the fast angle by performing a sequence of canonical transformations. In the following, this 
strategy will be translated in an explicit algorithm. Such a procedure will allow us to produce 
a final Hamiltonian satisfying two important properties: at the same time it provides a good 
approximation of the starting system, it only depens on the actions and one of the angles, i.e., 
it is integrable. 


2.2.1 Construction of the averaged normal form: the formal algorithm 

As discussed above, the normalization algorithm defines a sequence of Hamiltonians. This is 
done by an iterative procedure; let us describe the basic step which introduces starting 

from when both the values of the indexes ri and r 2 are positive. We assume that the 
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expansions of i) is such that 

ri — 1 


r2 —1 


s=o i>o 1=0 ( 8 ) 

l>r2 S>ri l>0 


where zj^'^ G Vi^sK V / > 0, 0 < s < ri, 


(n) 


G Vl,r,K V 0 < / < r2 , G Vpr,K 

Let us recall that the Hamiltonian ([5]) is suitable 


V / > r2 , ^ V Z > 0, s > ri 

to start the procedure with ri = r 2 = 1, after having set In formula (|^, one 

can distinguish the normal form terms from the perturbing part; the latter depends on p) in 
a generic way, while in the first two terms of ([5|) the fast variables can be replaced by the action 
r = + p ^)/2 . The (ri, r 2 )-th step of the algorithm formally defines the new Hamiltonian as 




(9) 


where exp ■ = X]j>o J\^x' series operator, with C^g = {g,x} (being {•, •} the 

classical Poisson bracket), g a generic function defined on the phase space and x generating 
function (for an introduction to canonical transformations by Lie series in the context of the 
Hamiltonian perturbation theory, see, e.g., IS]). The new generating function 
is determined so as to remove from the main perturbing terird //’’i/rj^its subpart that 
is not in normal form. This is done by solving the following homological equation with respect 

to xi7^ = xi7\p,^,Kv)- 


T (ri)Z, 

Xr 2 


( 0 ) 


I f(ri,r2-l-,ri) _ ^(ri) 

I Jr2 1 


( 10 ) 


where we require that is the new term in normal form, i.e. Zr7'^ = Zr7'^ (p, {7^ + p‘^)l2, A). 


Proposition 2.1 When Z^^ = (^^ + p‘^)j2 and ^77'''^ G 'Pr 2 ,r\K , then there exists a gen¬ 
erating function G Vr 2 ,riK ® normal form term Zr7'^ G Vr 2 ,riK solving the homological 
equation m- 


We limit ourselves to just sketch the procedure that can be followed so as to explicitly determine a 
solution of (|10p and, therefore, prove the statement above. First, we replace the fast coordinates 
(■^j v) with the pair of complex conjugate canonical variables {z, iz) such that ^ = (z —z)/\/2 and 
p = (z -\-'z)ly/2. Moreover, the homological equation (fT0|) has to be expanded in Taylor series 
with respect to {z,iz), using the slow coordinates (p, A) as fixed parameters (because they are 
not affected by the Poisson bracket L {rpZ^'^, since Z^^ do not depend on them). Therefore, 

Xr 2 

we solve term-by-term the equation (fTOll in the unknown coefficients Xm 2 ,m 3 {p, X) and Cmip,X) 
such that 

xi7Hp,z,X,iz)= [x^,,™3(p,A)z-^(m)-3] , z7^Hp,z,X,iz)=Y[Cmip,X)z"^7zr] ■ 

m2,m3 m 

^Let us recall that the size of p ‘‘£ Pi,bK is expected to decrease when the indexes s or r2 are 
increased, because the values of /r, p and s-re assumed to be small 
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At last, we express the expansions above by replacing {z,iz) for so as to obtain the final 

solutions in the form = Xr^ 2 ^\p^ p) and Zrl^'^ = (p, + p^)/2, A). 

The following property of the Poisson brackets is very useful for our purposes, and since the 
proof is immediate, it is omitted. 

Proposition 2.2 Let f and g be two generic functions such that f G Vr,sK o,nd g G 'Pr',s'K , 
then 

ii r + r' >2 {/, S'} G Vr+r',{s+s')K , else ^ {f,g} = 0 ■ 


In order to provide an algorithm easy to translate in a programming language, we are going to 
give explicit formulas for the new Hamiltonian and for its expansion that can be written 

as follows: 


ri — 1 


r2 


s=0 l >0 1=0 


I>r 2+1 


s>?’l l >0 


For the sake of simplicity in the calculation of we redefine the same quantity several 

times using the same symbol. This notation of the algorithm is more similar to its translation 
in a programming code, and, thus, more useful: let us introduce the recursive operation a 6, 
where the previously defined quantity a is redehned as a = a + 6. Therefore, we initially dehne 

j(ri,r 2 ,s) _ i,s) V / > r 2 when s = n or V / > 0 , s > ri. (12) 


Then, we consider the contribution of the terms generated by the Lie series applied to each 
function belonging to the normal form part as follows: 


f(ri,r2-,s+jri} J_^j ^(s) 
h+j{r2-2) ^ I 


y 1 <j < jf , 0</<//, 0<s<ri, 


(13) 


where the upper limits jj and If on the indexes j and I, respectively, are such that 

jf = 1 + 1 if r2 = 1 , jf = +00 if r2 > 2 , 

If = Too if s < ri , If = r2 if s = n , (14) 

k = r2 if s = ri , k = 0 if s > ri . 

For what concerns the contributions given by the perturbing terms making part of the expansion 
of FfTiG 2 -i) in ([8|), we have 


f{ri,r2-,s+jri) i_ ci f{ri,r2-Ls) 

Jl+j+ 2 - 2 ) ^ j! 


y 'i- <j < jf , l>k, s >ri , 


(15) 


where the limiting values for the indexes, that are jf and h , are dehned in (|14D . 

The redefinition rules (fTH|) and (fT^ are set so that the new perturbing part generated by the 
Lie series in ([9|) is coherently split in different terms according to their order of magnitude in p 
and their total polynomial degree in the actions. In fact, by applying repeatedly proposition 12.21 
to the redefinitions in (fT^ - (fT5]) . it is possible to inductively verify that ^ V I > 

h: s > ri. Therefore, the terms making part of the Hamiltonian in the expansion m 

share the same properties with those appearing in ([8]); this ensures that the normalization 
algorithm can be iterated so as to construct L/'(^i’'’ 2 +i)^ }firi,r 2 + 2 )^ ^ ^ ^ 
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2.2.2 Criteria for stopping the normalization algorithm, in order to perform a 
finite number of operations 

From an ideal point of view, we would be interested in producing the final Hamiltonian lim^^-s-oo bm^j-^oo 
where is defined as lim,. 2 _>.cxD V ri > 1. In fact, such a Hamiltonian would be 

integrable, because it would depend from the fast variables just through the action + ry^)/ 2 , 
due to the special form of the expansion (HID . In general, the problenH of the restrictions of 
domains prevents the convergence of the limits (with respect to both indexes ri and r 2 , when the 
standard sup-norm is used); thus, it is not possible to define the integrable Hamiltonian on any 
open set (see, e.g., 0 ). Actually, from a mathematical point of view, expansions of type m are 
asympotic series with respect to both ri and r 2 ; this means that we have to truncate the indexes, 
but in a way that optimize our result. For this purpose, we can proceed as follows. First, we 

introduce the functions Zbi>»' 2 ) = (^^ A) and 77.’^ 2 ) = 77bi'^2)(p, A, ry) that 

make explicit the splitting between the integrable and the perturbing parts in the expansion m, 
so that 


ri —1 




r2 




A»'i) 


Jl(ri,r2) _ 




s=0 l>0 1=0 l>r2+l s>ri l>0 

(16) 

Therefore, we look for the pair of upper indexes ( 7 ?i,i? 2 ) minimizing the sup-norm of 77 bi’'’ 2 ) 
(on the set of values of the variables that we are interested to study). This approach can be 
implemented with a suitable scheme of analytic estimates, in order to reduce exponentially the 
remainder with respect to the small parameters of the problem (see m, 0 and 0)- 

Such an optimal choice about the final values of the indexes ri and r 2 allows us to reformulate the 
algorithm, in such a way that it requires just R 1 R 2 normalization steps, constructing the finite se¬ 
quence of Hamiltonians ^ ^ ^ jp{Ri, 0 )^ ^ jj(R 2 ,Ri) ^ 

where FF^i+hO) = J 7 (n,R 2 ) V 1 < ri < F?i . 


While the algorithm has been rearranged so as to be performed in a finite number of normal¬ 
ization steps, it is evident that the redefinition rules reported in formulas (|13p and (USD would 
require to calculate infinitely many Poisson brackets. In order to avoid such a problem, we have 
to establish two truncation rules on the terms appearing in the expansions, so as to fix (a) their 
maximal exponent Smax related to the order of magnitude 0 {p^), (b) their total maximal degree 
Imax on the index I (that is equal to twice the degree in p plus the one in ^ and that in rj). If our 
formal algorithm is subject to a further restriction, such that it is limited to the calculation of 
functions belonging to classes of type VpsK with 0 < I < Imax and 0 < s < Smax , then it can be 
proved that it requires a finite total number of Poisson bracket^. Therefore, this newly restric¬ 
ted version of our algorithm is suitable to be translated in a programming code. In principle, 
the values of Imax and Smax should be chosen in order to optimize the final results; in practice, 
they are usually fixed (as well as the final indexes F 2 i and R 2 ) so as to fit with the available 
computational resources. 


®We emphasize that the most celebrated problem concerning the convergence of the normal forms, the accu¬ 
mulation of “small divisors”, does not affect our scheme, because the main integrable term in the homological 
equation cni, i.e. zf, depends just on the fast action -|- ?7^)/2. 

^Each of them needs a Hnite number of basic operations like derivatives, sums and products. 
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2.2.3 Approximate numerical integration based on the normalizing canonical trans¬ 
formation 


It is well known that Lie series induce canonical transformations in a Hamiltonian framework; 
this fundamental feature will allow us to design a numerical integration method, by using both 
the normal form discussed above and the corresponding canonical coordinates. In order to 
explicitly realize such a project, we have to introduce some more complicate notations. Let 
us denote with gg^ gf canonical coordinates related to the 

(?’i 5 ^ 2 )“th step. By appying the so called exchange theorem (see, e.g., [5]), we have that 

fj{ri,r2) ^p{ri,r2) ^(ri,r2) ^{ri,r2) p(ri,r2)'^ — jj{ri,r2-l) ^p(ri,r2) ^p(ri,r2) ^(ri,r2) ^{ri,r2) 

(17) 

where the variables related to the previous step, namely ’''2 i)^^hi ,»’2 i)^_\(n,r '2 i)^p{ri,r 2 
are given as 


(p(^l’^ 2 )(^(ri.r 2 )^^(ri.r 2 )^^(ri,r 2 )^^(ri,r 2 )^ = gxp f/l ^ (pirur 2 ) ^ ^irur 2 ) ^ ^{n,r 2 ) ^ piri,r 2 )\ . 

the r.h.s. of the equation above means that four Lie series must be applied separatedly to each 
variable, in order to properly define all the coordinates for the canonical transformation 
The whole normalization procedure can be described by the canonical transformation 


q(R2,Ri) ^ (^( 1 , 1 ) o 


0 (p(bR 2 )„ 


o ... O 


. . . o 


(19) 


Such a composition of all the intermediate changes of variables can be used for providing the 
following semi-analytical scheme to integrate the equations of motion: 


I 

piofi) (0), ^(0.0) (0), (0), (0)) 


(c(«i.« 2 ))’ 


piRi ,R2) (0)^ -R2) (0), (0), ’■^2) 

t 

Z{Ri,R2) 






CiRl,R2) 


^(Hl ,iJ2) ^ ,ii2) ^ ,K2) ^ ,fi2) 


( 20 ) 

where is the flow induced on the canonical coordinates by the generic Hamiltonian fC for 
an interval of time equal to t. Let us emphasize that the above integration scheme provides 
an approximate solution; from an ideal point of view (i.e., if all the expansions were performed 
without errors and truncations), formula (|20n would be exact if would correspond to 

the complete Hamiltonian On the other hand, is integrable and its flow is easy 

to comput^, reasons why using 2 ^(^i>^ 2 ) becomes valuable. The approximate solution provided 


®In order to explicitly describe the solutions of the equation of motions for the normal form Z^R'^^R^) ^ 
it is convenient to introduce the temporary action-angle variables such that _ 

V 2 r(^i’^ 2 ) cos£i^^’^^^ and where jg ^ constant of motion for the 

normal form (^p(Ri,R 2 ) ^ p(fli.fl 2 )^ ^(.Ri,R 2 )y gy considering ^g fixed parameter and 

using the standard quadrature method for conservative systems with 1 d.o.f., one can compute and 

\(Ri,R 2 )at any time t. The same can be done for the evolution of by evaluating the integral corres¬ 
ponding to the differential equation . For practical purposes, the application of the classical 

quadrature method can be replaced by any numerical integrator that is precise enough. Finally, the values of 
f(Ri,R 2 )u\ ggfi p(Ri,R 2 ) u\ ggg fig fifiectly calculated from those of the corresponding action-angle variables, that 
are r^R^’R2)(t) and 
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by the scheme (j 2 np is as more accurate as smaller the perturbing part jg with respect 

to (see their definitions in ()16p h 

As a hnal remark, we stress that it requires a finite total number of operations, if we limit the 
the expansions of the canonical transformations as discussed in the previous subsection. Thus, 
also the whole integration scheme () 20 p can be translated in a programming code. 

2.3 Tests on the accuracy of the normal form approximating the CPRTBP 
Hamiltonian 

In order to test the accuracy of our integrable normal form approximating the Hamiltonian, we 
have to choose a suitable surface of section for the comparison with the complete problem. The 
most logical choice about the sectioning of the flow is given by the surface defined by f = 0 , 
because f is a fast angle (recall the discussion about formula (|7P). 

In Fig. [T]we show the results of the comparison between the complete CPRTBP Hamiltonian 
and the normal form approximating it, in the case of the Earth-Moon system, which is dehned 
by its value of the mass parameter p = 0.01215058561. In the left panel, we show the surface of 
section numerically computed by considering the equations of motion related to the CPRTBP 
Hamiltonian. We take a set of 10 equispaced initial conditions, with p = 0, 4.188 < A < 4.45, 
i = 0. The value for P has been set in such a way that all the initial conditions keep the 
same value for the Jacobi constant as in L 4 . These orbits has been integrated up to recover 
1000 points over the surface defined by £ = 0. Those points are drawn in black in the space of 
variables (A,/)). 

For the same initial conditions as before, we have integrated the orbits also according to the 
semi-analytical scheme (f2n|) . up to collect also 1000 points over the surface. This computations 
have been automatically made by a code written in C, in such a way to use the expansions of 
the normal form and the canonical transformation which were preliminarly 

produced by using Mathematica. The truncations have been made according to the following 
values of the parameters ruling the extension of those expansions: Ri = 3, R 2 = 5, K = 5, Imax = 
5 and Smax = 3. These values imply reasonable computing times. The results are expressed in 
the middle panel of Fig. [TJ Finally, in the right panel, we show the comparison between the 
two surfaces of section. For all the orbits close to the equilibrium point L 4 , the correlation is 
extremely good, thought the approximation fails on reproducing exactly the tadpole orbits far 
from L 4 . 


3 Optimal transfers by means of integrable aproximations 

3.1 Baseline 

Since the new normal form has 1 d.o.f., it provides the integrable approximations of the motions 
of small bodies in the vicinity of the stable equilibrium points and suitable surfaces of section 
that are easy to compute. These tools can be used in many different studies. As explained in 
the Introduction, we focus in the design of maneuvres of small bodies (as spacecrafts or space 
debris) that are originally out of the stability region filled by tadpole orbits and we would like 
to situate into it. With such an aim, we design a method that allows to study the effectiveness 
of a set of impulses, done after considering a fixed initial position. 

As outline of the algorithm, we have to 
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Figure 1: Left panel, in black surface of section obtained by a numerical integration of the 
CPRTBP. Middle panel, the surface of section for the same initial conditions computed with the 
normal form approximating the CPRTBP Hamiltonian. Right panel, the comparison between 
the two surfaces of section. In the three plots, L 4 is located at (dvr/S, 0). 
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a. choose a starting point Pg, which corresponds to the position of the body to transfer at 
t = 0 and translate it to cartesian coordinates Pg = (xq, yo; 14 o) ^yo) • 

b. give an impulse 

= II AV|| cos /3 AI 4 , = II AV|| sin /? , (21) 

where ||AV|| and /3 are the modulus and direction of the impulse, respectively, in cartesian 
coordinates. This gives a new orbit, whose evolution in normalized variables can be com¬ 
puted with the integrable approximation of the Hamiltonian and checked in the surfaces 
of section. 

c. decide the acceptance or rejection of the impulse, according to whether it fulfills or not 
the conditions of the transfer. 

As a matter of fact, this mechanism can be applied for a physically adequate range of moduli 
and directions and not one orbit by one, in order to speed up the computations. 


3.2 Application and results 

We apply the method described above to the case of the Earth-Moon system. As starting 
point Pg , we choose a position slightly outside the real stability regiorU, estimated by numerical 
integrations of the full problem. Pg is given by A = 3.95, p = 0, ^ = 0.0885442 and rj = 0. For the 
transfer orbits, we apply 10“^ different impulses, for which 0 < ||AV|| < 0.01 and 0 < /3 < 27 r. 
Following §2.8 in i, m and (fT^ . we translate every orbit to normalized variables. In Fig. [2l 
we present the results of the translation of the new coordinates calculated after each impulse 
(for the first new point of the orbit on the surface of section) in the corresponding variables. 
The color scale represents the size of the impulse ||AV|| in the left panel and its direction /3 in 
the right panel. As expected, a stretching effect occurs while translating to the new system, and 
the impulses are more effective (i.e., same values imply a bigger difference with respect to the 
original point Pg) in some directions. 

In order to distinguish which impulses generate orbits that go deeper inside the stability 
region, we isolate the impulses pointing towards the inner part. We refine the values of the 
angles considered, taking 100 equidistant values /?, where 0.46 < /? < 0.86. We integrate their 
evolution with respect to the normal form and we find the curves described by these 

orbits on the surface of section. In the averaged normal form, each of these curves depicts the 
area associated to one of the actions defining the 2-D torus that is invariant with respect to the 
motion. Since as closer the torus is to the equilibrium point L 4 , as smaller this area is, we take 
as criterion for a good transfer a reduction of this value with respect to the initial one (Minimum 
Action criterion, since such an area corresponds to an action). In Fig. [3l we show the results 
of the computation of the area in two different color-scales. In the left panel, we can see for all 
the considered combinations of modulus ||AV|| and angles /?, the value of the area described 
by the transfered orbit. The lower border of the plot represents the orbits with ||AV|| = 0, so 
their areas (Aq) are equal and can be used as reference value. In the right panel, only orbits 
corresponding to areas < Aq are considered, the greater ones are fixed equal to Aq . This allows 
a better discrimination between all the orbits providing suitable impulses. 

®More specifically, Ps is in the chaotic region related to the stable/unstable manifolds emanating from L3. It 
has been fixed in order to be one of the points closest to L 4 among those belonging to the chaotic set and to the 
axis of the abscissas of the surface of section defined by ^ = 0. 
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Figure 2: Color scheme for the initial conditions produced after applying the impulses. Left 
panel, new positions in {X,p). The color scale represents the magnitude of the impulse ||AV||. 
Right panel, same as before with color scale representing the direction of the impulse /3. 



Figure 3: Left panel, values of the areas enclosed in the surface of section for every considered 
combination of ||AV|| and (3, the color scale represents the results of the calculations using the 
averaged normal form. Right panel, same as before where the computed areas are reported only 
if they are smaller than the value (^o) corresponding to the original initial condition Pg . 
















Maneuvers based on integrable approximation 
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Both plots shows that the selection of the angle for the impulse is not trivial, since small 
differences generate very different results. Consider, for instance, the impulses with (3 = 0.50 
and with (3 = 0.55, for ||AV|| > 0.006). The former give the best approach to the optimal 
transfer we look for, while the latter generate orbits which are actually further from L 4 than 
that related to the initial condition Pg . In general, since we take small values for the impulses, 
in the correct directions, bigger ||AV|| implies smaller final areas. The best choices for the 
transfer correspond to the impulses for which (3 € [0.50, 0.68] , and suitable values for their 
sizes, following the darker stripes in Fig. [3j 


4 Conclusions 

In this work we explicitly construct an integrable normal form approximating the CPRTBP 
Hamiltonian. This is done by reformulating the approach described in [3] so as to use some 
more modern techniques, developed in the framework of the Hamiltonian perturbation theory 
and mainly based on the Lie series formalism. This allow us to design an algorithm that can be 
fully translated in programming codes. In particular, we produce a truncation of the normal form 
whose the expansion in (II 6 D highlights that it is an average of the CPRTBP Hamiltonian 
with respect to the angle associated to its fast dynamics. The first results provided by this 
revisited approach are encouraging: in some suitable surface of sections our algorithm provides 
very good approximations of the tadpole orbits close enough to L 4 — L 5 . However, we are aware 
of the fact that the accuracy of our expansions must be strongly improved in order to face 
challenging concrete problems in a region far from those equilibrium points. In our opinion, the 
main constraint on the quality of our results is due to the truncations on the Fourier series in 
the slow angle A. We think that this limitation can be removed, by representing the dependence 
on A in a suitable way, so as to avoid Fourier expansions. We plan to investigate such a new 
approach in the near future. 

Furthermore, we show a first astrodynamical application starting from our calculation of 
the integrable normal form which approximates the CPRTBP Hamiltonian. We design 

an algorithm that allowed to compute optimal transfers between orbits in the neighborhood 
of the equilateral Lagrangian equilibrium points. We generate impulses in cartesian variables 
according to formula (I2ip . on a grid of values for ||AV||, the magnitude of the impulse, and (3, 
the direction. For those, we are able to discriminate the suitable transfers, that imply a final 
orbit closer to the equilibrium point, according to a new criterion. Using our normal form as an 
approximation of the complete CPRTBP, we estimate the area enclosed by the final orbit, and 
minimizing this quantity, we select the best candidates for the transfer. A careful inspection of 
the plots shows that the best candidates are highly depending on the size of the impulse, and 
very sensitive to the changes on the angle /3. 
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